// QUIRKSIM
// Copyright (c) 2009 Gregory Andrew Robinson (garobinson@ucdavis.edu)
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.

#ifndef MAIN_CC_
#define MAIN_CC_

#include <stdexcept>
#include <gflags/gflags.h>

#include "quirksim.h"

DEFINE_double(dt, 0.01, "Time step for each iteration.");
DEFINE_double(max_time, 1, "Time to spend during simulation.");
DEFINE_int32(num_points, 100, "Number of string parameterization points.");
DEFINE_int32(plot_every, 1, "Number of iterations for each log event.");
DEFINE_double(x_min, -10, "X comp of first quirk position.");
DEFINE_double(x_max,  10, "X comp of second quirk position.");
DEFINE_double(sigma, 10, "String \"tension\".");
DEFINE_double(point_separation, 0.1, "String \"tension\".");
DEFINE_double(height, 1, "Height of initial waveform.");

double triangle(double x, double x_min, double x_max) {
  return 1.0 - fabs(x) / (x_max - x_min) * 2.0;
}

double sine(double x, double x_min, double x_max) {
  return sin(x * 2 * PI / (x_max - x_min));
}

int main(int argc, char *argv[]) {
  google::ParseCommandLineFlags(&argc, &argv, true);

  double dt = FLAGS_dt;
  double max_time = FLAGS_max_time;

  unsigned int num_points = FLAGS_num_points;
  unsigned int plot_every = FLAGS_plot_every;

  double x_min = FLAGS_x_min;
  double x_max = FLAGS_x_max;

  double sigma_ = FLAGS_sigma;
  double point_separation_ = FLAGS_point_separation;

  QuirkSystem sys(sigma_, point_separation_);

  sys.mutable_points()->push_back(new FixedEnd(Hep3Vector(x_min, 0, 0)));

  double x;
  for (unsigned int i = 1; i <= num_points; i++) {
    x = x_min + i * (x_max - x_min) / (num_points + 1);
    sys.mutable_points()->push_back(
        new StringPoint(
            Hep3Vector(x, FLAGS_height * sine(x, x_min, x_max), 0),
            Hep3Vector(0, 0, 0),
            sys));
  }

  sys.mutable_points()->push_back(new FixedEnd(Hep3Vector(x_max, 0, 0)));
//  sys.mutable_points()->push_back(
//    new Quirk(Hep3Vector(x_max, 0, 0),
//              Hep3Vector(),
//              0.000001,
//              sys,
//              true));

  int i = 0;
  for (std::vector<ParamPoint*>::const_iterator it = sys.points().begin();
       it < sys.points().end();
       it++) {
    if (i % plot_every == 0) {
      Print((*it)->x());
      printf("\n");
    }
    i++;
  }

  for (double elapsed = dt; elapsed <= max_time; elapsed += dt) {

    sys.RK4Step(dt);
  }

  i = 0;
  for (std::vector<ParamPoint*>::const_iterator it = sys.points().begin();
       it < sys.points().end();
       it++) {
    if (i % plot_every == 0) {
      Print((*it)->x());
      printf("\n");
    }
    ++i;
  }
}

#endif  // MAIN_CC_
